Dissipative flow and vortex shedding in the Painleve boundary layer of a Bose 
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Raman et al.[l] have found experimental evidence for a critical velocity under which there is no 
dissipation when a detuned laser beam is moved in a Bose-Einstein condensate. We analyze the 
origin of this critical velocity in the low density region close to the boundary layer of the cloud. In 
the frame of the laser beam, we do a blow up on this low density region which can be described by a 
Painleve equation and write the approximate equation satisfied by the wave function in this region. 
We find that there is always a drag around the laser beam. Though the beam passes through the 
surface of the cloud and the sound velocity is small in the Painleve boundary layer, the shedding 
of vortices starts only when a threshold velocity is reached. This critical velocity is lower than the 
critical velocity computed for the corresponding 2D problem at the center of the cloud. At low 
velocity, there is a stationary solution without vortex and the drag is small. At the onset of vortex 
shedding, that is above the critical velocity, there is a drastic increase in drag. 

PACS numbers: 03.75.FL02.70.-c 



Dilute Bose-Einstein condensates have recently been 
achieved in confined alkali-metal gases and the study 
of vortices therein is one of the key issues. Raman et 
al. [1, 3], Onofrio et al. [2] have studied dissipation 
in a Bose Einstein condensate by moving a blue detuned 
laser beam through the condensate at different velocities. 
They found experimentally a critical velocity for the on- 
set of dissipation. This critical velocity has been related 
to the one found by Frisch et al. [4] for the problem of a 
2D superfluid flow around an obstacle in the framework 
of Nonlinear Schrodinger Equation (NLS): below a crit- 
ical velocity, the flow is stationary and dissipationless, 
while beyond this critical velocity, the flow around the 
disc becomes time dependent and vortices are emitted. 
Numerical simulations have been done for this type of 
problem in 2D [5] and 3D [6, 7]. In particular, the direct 
3D simulation of [7] shows the plot of the drag against the 
velocity. A critical velocity can be numerically computed 
when the drag becomes nonzero, but no precise mecha- 
nism of vortex nucleation is described by the authors. 
This critical velocity has been analyzed theoretically for 
a homogeneous 2D system [8] and an inhomogeneous 2D 
system [9, 10]. 

In this paper, we want to take into account the 3D ge- 
ometry of the experiment of [1-3]. Our aim is to under- 
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stand the mechanism of vortex nucleation in the bound- 
ary region. Indeed the analysis of [4] allows to understand 
what is happening in the interior of the cloud, where the 
kinetic energy is negligible in front of the interaction en- 
ergy. In the region where the laser beam crosses the 
boundary of the cloud, the sound velocity gets small, 
since the amplitude of the wave function becomes small. 
There, the kinetic energy term can no longer be neglected 
in front of the trapping and interaction terms. We blow 
up this region in such a way that the trapping potential 
varies linearly with the distance to the boundary and far 
away from the laser beam, the wave function is then given 
by a Painleve equation. We analyze the behavior of the 
wave function in the frame of the laser beam. The real 
experiments are quite complex, and in particular here we 
do not take into account the oscillations and acceleration 
of the beam but we believe that our analysis allows to 
understand the mechanism of increase of drag. One of 
our main results is that there is always a drag around 
the laser beam and this drag grows continuously. At low 
velocity, the drag is not a consequence of the shedding 
of vortices, and finally of a time dependent density and 
velocity field. The origin of this drag is in the radiation 
condition for the wavefield: the motion changes contin- 
uously the structure of the solution seen in the frame of 
reference of the "fluid" at infinity. We study the transi- 
tion toward a time dependent regime of vortex shedding, 
which happens at a critical velocity. The critical velocity 
that we find is lower than the 2D critical velocity at the 
center of the cloud coming from the computation of [4]. 
Vortices are nucleated close to the boundary of the cloud 
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and the tubes grow and detach to form rings that move 
downstream. When tubes are emitted, significantly large 
drag values are observed. The drag increases smoothly 
as the velocity increases. 

The dynamics can be modeled using the Gross 
Pitaevskii equation at zero temperature with an exter- 
nal trapping potential V tr = m/2{uj 2 x 2 + uj 2 y 2 + u 2 z 2 ). 



h 

ihdf^ = A^p 

2m 



(V tr + Ng\*\ 2 )*. 



If an object is moved inside the condensate, V tr has to 
be replaced by V tr + V b, where V Q b depends on x — vt. 
Based on the experimental data of [1, 2], we take a = 
mg/Anh 2 = 2.94 ram, N = 1.2. 10 7 , u y = uj z = 377s- 1 , 
and uo x = \uo Zl with A = 0.3. We also define the char- 
acteristic length d = (h/mcUz) 1 / 2 = 2.71 pm and a small 
nondimensionnalized parameter e given by 
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We find that e = 6.21 10 -3 which may be viewed as 
small parameter and allow rescaling the equation near 
the edge of the condensate. Re-scaling the distances by 
R = d/y/e = 34.4 /ira, the time by l/(eu z ), we have 
^(r, t) = R 3 / 2 (j)(r , t) where f = Rr. In these new units, 
the radii of the condensate are R y = R z = 0.65 and R x — 
2.18. The laser beam is modeled by an obstacle which 
is a cylinder C of axis z and radius I = 0.19 on which 
ip = 0. We will work in the frame where the obstacle 
is stationary. Outside the obstacle, the equation can be 
rewritten as 

-2id t ^ = A^ + ^(ptf-M 2 )^ 

where ptf = Po — (A 2 x 2 + y 2 + z 2 ) is the Thomas 
Fermi limit density and po — 0.42 is the rescaled chem- 
ical potential. Note that |^| 2 is close to its Thomas 
Fermi value ptf except near the obstacle and near the 
boundary of the cloud. This boundary layer has a thick- 
ness of order e 2 ! 3 so that we rescale the domain with 
tp(x,y,z) = e 1 / 3 u(x, y, z), where x = xje 2 ! 3 ^ y = yje 2 ! 3 
and z = (yfpo — z)/s 2 / 3 , v = ve 2 ! 3 . By blowing up 
the boundary of the cloud near z = 0, and truncating 
at z = L, the rescaled layer thickness, we see that the 
modulus of the stationary solution in the boundary layer 
for \x\ and \y\ large, that is far away from the obsta- 
cle, is given by the solution of the first Painleve equation 
[11, 12]. 

p" + (2z^-p 2 )p = 0, p(0) = 0, p(L) = ^2y^L. (1) 



We choose the size of the boundary layer L so that 
e 2 l 3 L = 3^/po/lO. This is based on the consideration 
that, on the one hand, e 2 ^ 3 L should be suitably small so 
that 2zy/po is a good approximation for ptf — Po — z 2 
in the boundary layer and on the other hand the critical 
velocity at z = L is not too different from the critical 



velocity at the center of the cloud. The obstacle is now 
a cylinder of radius a = Ije 2 ! 3 = 5.6. 

The obstacle moves at the rescaled velocity v = 
/(e 1 / 3 ^^), and in the frame of the obstacle, the 



u exp 

equation becomes 

-2id t u = Au - 2ivd x u + (2z^ - \u\ 2 )u. (2) 

We want to understand the behaviour of solutions de- 
pending on v. If we restrict (2) to z — L, we can perform 
a similar analysis to [4] and get the value of the crit- 
ical velocity for the onset of vortex shedding and find 
v 2 = 2y / poL/ll = 2c 2 /ll, where c s is the sound velocity. 
Of course, we cannot apply this analysis in the low den- 
sity region, since there the sound velocity gets close to 0. 
Another mechanism has to be understood. The rescaled 
drag around the obstacle is 



drag = ^ J (u x u n - u x u n ) dl dz. 



(3) 



We first analyze the stationary solution of (2) in the 
very low density region, where the system is very dilute 
and one can neglect the nonlinear term. In fact, a pre- 
cise condition is that p 2 is less than v 2 , which gives a 
truncation point z c at which p 2 (z c ) = v 2 . It is rather 
straightforward in classical scattering theory to compute 
the perturbed wavefield and finally the drag on the ob- 
stacle (a related problem, the scattering of sound by a 
cylinder, is treated in [13]). In the low density region, it 
is reasonable to look for u with the following ansatz 



u(x,y,z) =p(z)ilj(x,y)e v 



(4) 



We can first approximate p(z) in this region by an Airy 
function given by the solution of p" + 2zp^/po = 0, that 
is, by defining z 3 = l/(2^/joo), we have 

P(z) « V2Ai(^) « -I=(Z£)l e ^(-|(^)l). ( 5 ) 
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Then, outside the obstacle, ip is a solution of the 2D 
Helmholtz equation 



Aip + v 2 ^) = 



(6) 



with ip = for r = a, the obstacle boundary, and ip w 
e lvx at infinity. This solution can be computed [13] in 
terms of Bessel functions Jk and 7V fc . One finds that, at 
leading order for v small, the 2D drag of ip is proportional 
to 

v 2 JS(v)J 1 (v)N 2 (v)/NS(v) - v/\og 2 v. (7) 

The total drag has to be multiplied by the integral of p 2 
along the z axis to the truncation point z c defined by 
p 2 {z c ) = v 2 . Direct calculation gives 



v f Zc , 
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2 dz^C^ r „ 



(8) 



In conclusion, the total scattering drag tends to zero at 
low speed. It is plotted in Figure 1 (solid line). 

To understand how the solutions of (2) behave, we nu- 
merically integrate the equation in a computational do- 
main of dimension 60 x 60 x L with periodic boundary 
conditions in x and y and taking u = on the boundary 
of the obstacle and away from the condensate (z = 0). 
At the truncated surface z = L inside the condensate, 
we use the condition (d / dz)(u/p) = 0, where p is the 
solution to the Painleve equation (1) described before. 
The numerical solution is computed based on a continu- 
ous piecewise quadratic finite element approximation in 
space and the Runge-Kutta fourth-order in time integra- 
tion scheme. Using p as initial condition, we first com- 
pute the solution of (2) for some time by adding a damp- 
ing coefficient, that is, we replace iut in (2) by iut(l + ij). 
For small velocity, this effectively drives the numerical 
solution of (2) close to a stationary solution. Then, we 
continue the integration with a much reduced damping 
coefficient 7 = 0.02 or with no damping 7 = 0. 

In what follows, we will divide the velocity by the 
sound velocity at center c s = ^/2po/e 1 ^ 3 . In Figure 1, 
we plot the drag vs. the velocity divided by c s . For 
a given velocity value, the drag is the value obtained 
through time averaging of (3). We have verified that with 
different small values of 7, the drag calculation remains 
essentially the same. 
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FIG. 1: Drag vs. v/c s : for (8), — o— for numerical solu- 
tion of (2); insert: zoomed in for small v. 

For v small, we find that the solution is almost sta- 
tionary. Surface oscillations are present near z = 0, and 
the drag is small, but not zero. See Figure 2 for plots of 
the solution. The drag computed in this regime fits very 
well to the cubic growth given by (8). 

When v is increased, at a critical velocity v c /c s « 0.2, 
the surface oscillations develop into small handles that 
move up and down the obstacle without detaching; see 
Figure 3. 

There is no stationary solution, but no vortex shedding 
either: the small handles move up the obstacle to a crit- 
ical z value and down. This instability may be related 
to the one discussed by Anglin [14]: in our scaling, the 
critical velocity found in [14] is 0.2. At this stage, the 
solutions do not produce large drag nor vortex shedding. 

It is only for larger velocities (v/c s > 0.25) that the 
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FIG. 2: Isosurface snapshot of \u\: surface waves for v = 0.08 
and v — 0.2. 




FIG. 3: Isosurface snapshots of \u\ at different times for v = 
0.24: formation of vortex handles. 



handles move up to the top, detach from the obstacle 
and produce significant drag. This is a wholly nonlinear 
phenomenon and most likely cannot be described by a 
linear analysis. 

Let us describe the solutions for v/c s > 0.25 illustrated 
in Figure 4. The vortex handles seem to first nucleate 




FIG. 4: A sequence of isosurface snapshots of \u\ for v = 0.28: 
a) formation of vortex handles, b) detachment from obstacle, 
c) bending of vortex tubes and d) formation of vortex half 
rings. 

near z = and are top connected to the obstacle. As 
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time increases, the bottom ends move away from the ob- 
stacle in a slightly down stream direction while the top 
end moves up along the obstacle (Figure 4a). When the 
top ends of the vortices become close to z — L, the bot- 
tom ends reverse their trend of moving away from obsta- 
cle. Instead, they move back to the bottom of the ob- 
stacle, as if the handles prefer certain curvature (Figure 
4b). Eventually, the top ends of the handle move away 
from the obstacle and produce a pair of vortex tubes 
with their bottom ends at the bottom of the obstacle 
(Figure 4c). The handles merge into a half vortex ring, 
this half ring moves both upward and downstream (Fig- 
ure 4d). Near z = 0, the solution can be approximated 
by the solution (4) and this solution does not have vor- 
tices, so the instability creates the vortex but the vortex 
moves away. Vortex detachment happens only at suffi- 
ciently high density, in the region where the nonlinear 
term in the equation dominates. The direction of the 
vortex displacement is due to the velocity of the flow and 
the self interaction of the vortex on itself, which gives a 
movement along its normal vector. Meanwhile, while the 
vortex ring starts to detach from the obstacle, another 
pair of vortex handles is forming near the obstacle. The 
above process repeats itself. Note that we have truncated 
the domain close to the boundary of the cloud, so that 
the half ring we compute would correspond to a closed 
ring in the experiments. 

We have to point out that the critical velocity we have 
found for the onset of vortex shedding is lower than the 
critical velocity for the 2D problem at z = L. In this 
case V2d/c s = 0.35. So the inhomogeneity in the con- 
densate lowers the critical velocity from the 2D value. 
One can check that for different L, the critical velocity 
does not change. This is verified by our numerical com- 
putation where we have used two boxes with one about 
50% higher in z than the other, and there is little change 
in the drag plots, nor there is any significant difference 
in the dynamic behavior of the solutions. 



In the experiments [1-3], the drag is plotted vs velocity 
and a critical velocity can be defined when a sharp slope 
is observed in the drag plot. The critical velocity in [1] 
is very similar to ours, though slightly smaller. This is 
certainly due to the finite extent of the condensate in the 
x, y direction. Indeed, our simulations have not taken 
into account that the cloud is narrower in the y direction 
than along the x. We can check that for the 2D problem, 
this geometry lowers the velocity. On the other hand, 
our computations indicate that the inhomogeneity in the 
z direction and the soft boundary of the laser beam are 
well accounted for by our problem. 

Summary. We have studied the onset of dissipation 
in the Painleve boundary layer of a BEC when a detuned 
laser beam is moved in the condensate. We do a change 
of frame and blow up the low density region near the 
boundary of the cloud to write the equation for the wave 
function in this region: z = is now the boundary of 
the cloud and z large is the center. For small velocity, 
there is a drag around the obstacle due to radiation, but 
no vortex is generated: it is a stationary flow, which is 
supersonic near z = 0, but subsonic for z larger. On the 
other hand, when the critical velocity is reached, the in- 
stability propagates towards the top, a vortex handle is 
nucleated and detaches from the obstacle to form vortex 
rings that move away. Our aim was to understand the 
origin of vortex shedding. The critical velocity is lower 
than for the 2D problem. There is a drag for all veloc- 
ity, it increases smoothly with the velocity, and there a 
significant increase at the onset of vortex shedding. 
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